### Replication script for Tables 1 and 2 ###
library(tidyverse) # 1.3.2
library(fixest) # 0.10.4

# Load in
load("lb_clean.RData")
load("district_covars.RData")

# Table 1
# Combine survey data with district covariates
lb_all <- left_join(lb_clean, district_covars)

presence_only <- lb_all %>%
  filter(crim_presence == 1)

m1 <- feols(crim_presence ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_all)
m2 <- feols(crim_presence ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_all)
m3 <- feols(crim_presence ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = lb_all)
m4 <- feols(crim_gov ~ confidence_government + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only)
m5 <- feols(crim_gov ~ local_gov_responsive + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only)
m6 <- feols(crim_gov ~ police_corrupt + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs | district_code, weights = ~ wt, data = presence_only)

etable(m1, m2, m3, m4, m5, m6,
       order = c("Confident in government", "Local gov. is responsive", "Police are corrupt"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(confidence_government = "Confident in government",
                local_gov_responsive = "Local gov. is responsive",
                police_corrupt = "Police are corrupt",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs"),
       tex = T)

# Table 2
m1 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only)
m2 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = presence_only)
m3 <- feols(crim_gov ~ coercive_presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_all)
m4 <- feols(crim_gov ~ presence_pca + SEX + AGE + EDATTAIN + UNEMPLOYED + salary_covers_needs + road_density + log_population + pop_density + luminosity_pc | idenpa, weights = ~ wt, data = lb_all)

etable(m1, m2, m3, m4,
       order = c("State presence \\(coercive\\)", "State presence \\(all\\)"),
       digits = "r3",
       cluster = ~ district_code,
       digits.stats = 2,
       signif.code = c("**"=0.01, "*"=0.05, "."=0.10),
       dict = c(coercive_presence_pca = "State presence (coercive)",
                presence_pca = "State presence (all)",
                SEXMale = "Male",
                AGE = "Age",
                EDATTAIN = "Educational attainment",
                UNEMPLOYED = "Unemployed",
                salary_covers_needs = "Salary covers needs",
                road_density = "Road density",
                log_population = "ln(Population)",
                pop_density = "Population density",
                luminosity_pc = "Luminosity per capita"),
       tex = T)


